if !d.name eq 'PS' then begin
    device,xsize=18,ysize=12,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
;!p.multi=[0,3,2]
;!p.charsize=3.5
!p.multi=[0,3,1] & !p.charsize=2.0 & !x.margin=[3.75,0] & !y.margin=[3,0.5]

restore,filename='./khelicity.sav'
restore,filename='./mhelicity.sav'

thg=where(th GE -89 AND th LT 89)
rgood=0.95
dd=where(r le rgood)
nlev=64

xlength=0.26
xinterval=0.06
xdelta=0.04
xpos1=xinterval
xpos2=xinterval+xlength

; kinetic helicity
tkhel=transpose(khel)
lev=grange(min(tkhel[dd,*]),max(tkhel[dd,*]),nlev)
print,'levels R_\phi\theta = ',minmax(lev)
cpos=[xpos1,0.15,xpos2,0.96]
bpos=[xpos1+xdelta,0.065,xpos2-xdelta,0.08]
contour,smooth(tkhel[dd,*],1),x[dd,*],y[dd,*],/fi,nl=64,/iso,levels=lev, $
	title='!8h!Dk!N=<u`.!4x!8`>!6',pos=cpos,charsize=2,ystyle=4,xstyle=4
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.72),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=3,li=0
oplot,[0,0],[0.61,0.96],thick=3,li=0
colorbar,range=[min(lev),max(lev)],pos=bpos,/horizontal,xtickformat='(E9.1)'$
        ,xticks=2,xtickv=[min(lev),0.,max(lev)],yaxis=0,char=1.7
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
; kinetic helicity
tmhel=transpose(mhel)
lev=grange(min(tmhel[dd,*]),max(tmhel[dd,*]),nlev)
print,'levels R_\phi\r flux = ',minmax(lev)
cpos=[xpos1,0.15,xpos2,0.96]
bpos=[xpos1+xdelta,0.065,xpos2-xdelta,0.08]
contour,smooth(tmhel[dd,*],1),x[dd,*],y[dd,*],/fi,nl=64,/iso,levels=lev,$
	title='!8h!Dm!N=<b`.a`>/!4ql!8!D0!N!6',pos=cpos,charsize=2,ystyle=4,xstyle=4
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.72),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=2,li=0
oplot,[0,0],[0.61,0.96],thick=2,li=0
colorbar,range=[min(lev),max(lev)],/horizontal,pos=bpos, $
         xtickformat='(E9.1)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.7
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
; total alpha effect
alpha=-tkhel+tmhel
lev = 2.*max(abs(alpha))*indgen(nlev)/float(nlev-1)-max(abs(alpha))
lev = grange(min(alpha[dd,*]),max(alpha[dd,*]),nlev)
print,'Total R_\theta\r = ',minmax(lev)
cpos=[xpos1,0.15,xpos2,0.96]
bpos=[xpos1+xdelta,0.065,xpos2-xdelta,0.08]
contour,smooth(alpha[dd,*],1),x[dd,*],y[dd,*],/fi,nl=64,/iso,levels=lev,$
	title='3!4a!8/!4s!8=-h!Dk!N+h!Dm!N!6' ,pos=cpos,charsize=2,ystyle=4,xstyle=4
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.71),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=2,li=0
oplot,[0,0],[0.61,0.96],thick=2,li=0
colorbar,range=[min(lev),max(lev)],pos=bpos,/horizontal,$
         xtickformat='(E9.1)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.7
;
;partvelvec,fx[*,*],fy[*,*],x[*,*],y[*,*],/over,fraction=0.1,length=0.08

!p.multi=0
end
